SANN: A parameter-free nearest-neighbor algorithm 
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We propose a parameter-free algorithm for the identification of nearest neighbors. The algorithm 
is very easy to use and has a number of advantages over existing algorithms to identify nearest- 
neighbors. This solid-angle based nearest-neighbor algorithm (SANN) attributes to each possible 
neighbor a solid angle and determines the cutoff radius by the requirement that the sum of the 
solid angles is 4ir. The algorithm can be used to analyse 3D images, both from experiments as well 
as theory, and as the algorithm has a low computational cost, it can also be used "on the fly" in 
simulations. In this paper, we describe the SANN algorithm, discuss its properties, and compare it 
to both a fixed-distance cutoff algorithm and to a Voronoi construction by analyzing its behavior 
in bulk phases of systems of carbon atoms, Lennard- Jones particles and hard spheres as well as in 
Lennard- Jones systems with liquid-crystal and liquid-vapor interfaces. 



I. INTRODUCTION 

In most studies of many-particle systems, one is con- 
fronted with the task of determining the nearest neigh- 
bors of a particle, or set of particles. Interestingly, while 
identifying nearest neighbors is an important component 
of various analyses, and is sometimes even needed to eval- 
uate intereraction potentials, there is no unique definition 
of a nearest neighbor and, as a result, what one defines 
as a nearest neighbor is typically dependent on the ques- 
tion at hand. The two most common algorithms for de- 
termining nearest neighbors are i) a fixed-distance cutoff 
and ii) a Voronoi construction However, many exten- 
sions, and other definitions have been used as well, see 
for instance Refs. [U, HI- 

A fixed-distance cutoff is the obvious choice in simu- 
lations with particles interacting through a short-range 
potential, where each nearest neighbor is an interaction 
partner and the cutoff distance corresponds to the in- 
teraction range. Additionally, fixed-distance cutoffs have 
also been used in determining nearest neighbors for struc- 
tural analyses such as calculatingbond-order parameters 
in nuclcation studies (e.g. see [4[). However, in these 
cases, the "fixed-distance" is not well defined. Arguably, 
the first minimum of the pair correlation function g(r) 
(also known as radial distribution function) is a reason- 
able choice for the cutoff, as it relates to the neighbors 
in the first coordination shell. However, the precise loca- 



tion of this minimum depends on both the system's de- 
tails and thermodynamic conditions and therefore must 
be determined every time either one is changed. Addi- 
tionally, the cutoff is defined for the entire system and, 
as such, is not appropriate for systems with large density 
gradients, such as occur naturally in nucleation studies, 
in systems in the presence of gravity or in systems with 
interfaces. 

In contrast, a Voronoi construction [l| is based on 
purely geometric constraints and is parameter free. In 
addition to identifying nearest neighbors, this method 
can be used to determine geometric properties like edges 
and faces shared between these neighbors - data that are 
frequently useful for structural analysis and classification. 
Based on the local environment around a particle, it is 
more appropriate than a fixed-distance cutoff in the case 
of density gradients. However, there are also a number of 
inherent problems with a Voronoi construction, some of 
which will be discussed in this manuscript. First of all, 
the method is computationally expensive and hence is 
rarely used on-the-fly in simulations. More importantly, 
it is not robust against thermal fluctuations. In a crys- 
tal, thermal fluctuations which cause particles to fluctu- 
ate around their equilibrium lattice sites can spuriously 
increase the number of particles which share a small face 
with the target particle [5, 6] and hence increase the num- 
ber of particles identified as nearest neighbors. There ex- 
ist extensions to the Voronoi construction which aim to 



2 



increase the robustness against these fluctuations 0- 
ITol |. however, they typically introduce parameters, re- 
moving the "parameter free" advantage of the algorithm, 
and they further increase the computational cost. 

Looking at the advantages and disadvantages of both 
the fixed-distance cutoff and the Voronoi algorithm, we 
suggest a list of features which a "good" nearest neighbor 
algorithm should have. Specifically, an algorithm should 
i) be able to deal with systems with inhomogeneous den- 
sity, ii) be stable against thermal fluctuations, iii) be pa- 
rameter free and iv) be computationally inexpensive. In 
this manuscript, we propose, with these goals in mind, a 
simple algorithm for the identification of nearest neigh- 
bors: the solid-angle based nearest-neighbor algorithm 
(SANN). This method is based on similar principles to 
a theory used by Corwin et al. [IlT ] which used solid 
angles to predict the number of nearest neighbors. It 
is also similar to a Voronoi construction as it does not 
require tuneable parameters. However, SANN is com- 
putationally significantly less expensive than a Voronoi 
construction. In fact, its computational cost only slightly 
exceeds that of a fixed-distance cutoff making it suitable 
for on-the-fly use in simulations. In order to compare 
our algorithm with the fixed-distance cutoff and Voronoi 
construction, we apply all three methods to monodisperse 
hard spheres, Lennard- Jones liquid and fee crystal bulk 
phases, the 3-fold coordinated liquid carbon and graphite 
phases and the 4-fold coordinated liquid carbon and di- 
amond phases. Additionally, to examine the behavior of 
the SANN method at interfaces we also analyze a liquid- 
crystal and a liquid- vapor Lennard- Jones two-phase sys- 
tems. 



II. METHOD 
A. Description of SANN method 

As mentioned in the introduction, there exists no 
unique definition of a nearest neighbor. Consequently, 
it comes as no surprise that our SANN algorithm intro- 
duces a definition that differs from those of existing algo- 
rithms. Yet it has many similarities with the definitions 
of the fixed-distance cutoff and the Voronoi construction, 
as it is based on similar concepts. 

Consider a dense system with excluded volume, where 
we have a particle i located at position r,; surrounded 
by particles {j}. The fixed-distance cutoff defines the 
nearest neighbors of particle i to be all the particles of 
{j} with a distance to i smaller than the cutoff-distance. 
However, as mentioned in the introduction, the problem 
with this definition is in choosing that distance. This is 
where our SANN algorithm comes into play. For each 
particle i SANN determines an individual cutoff distance 
which we call the shell radius. It depends on the 
local environment of particle i and includes its m nearest 
neighbors. Since the cutoff distance is now a local prop- 
erty, the algorithm is suitable for systems with inhomoge- 




Figure 1: Definition of the angle 6ij associated with a neigh- 
bor j of particle i. Here, nj is the distance between both 



particles and i?' m ' is the neighbor shell radius. 



neous densities. For the computation of r\ SANN uses 
a purely geometrical construction, as does the Voronoi 
tessellation. Thus, the algorithm is parameter-free and 
scale-free. In the following we describe the geometrical 
construction and how m and R^ are determined. 

First, we assume the particles {j} surrounding i are 
known and ordered such that n t j < for all j. This 

relates R^™ 1 and m in the following manner: 



r- < 7? (m) < r _li 



(i) 



Then, starting with the particle closest to i we associate 
with each potential neighbor j an angle 6ij based on the 
distance between the particles rj j = \ fj — r\\ and the yet 



undetermined shell radius R^ as depicted in Figure [TJ 
SANN defines the neighborhood of a particle i to con- 
sist of the nearest (i.e. closest) m particles {j} such that 
the sum of their solid angles associated with 6i t j equals 
4ir, i.e. 



47T 



2tt[1 



cos 



3=1 



2^(1-^-/^). (2) 



We point out that while the number m and the shell 
radius i?j are not known yet, they are not independent: 
once one is known it is straightforward to determine the 
other. Also note that since the solid angle contribution 
for a single neighbor is always less than 2-7r, m must be 
at least 3. 

To visualize this idea imagine each solid angle as a cone 
with its apex point located at particle i and the cone's 
base center located at neighbor j. For a complete set of 
nearest neighbors all those cones stack to fill a spherical 
volume around i with a radius corresponding to the shell 
radius iij • Obviously, cones are not space-filling (i.e. 
they don't stack without gaps), so for the sum of solid 
angles to equal Att some cone overlap does occur. 
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Figure 2: 2d comparison of the Voronoi construction (left panel) and the SANN algorithm (right panel). In both panels the 
nearest neighbors are indicated with red lines. To facilitate the comparison between a Voronoi construction and the SANN 
algorithm, we make use of the fact that the algorithm is scale free, i.e. ri^/R^ = (0.57 , i I j)/(0.51£j ) where 0.5rij is simply 
the distance from particle i to the midpoint between i and j, and 0.5i?' m ' is half the shell radius. The green circles have a radius 
equal to the half the shell radius of the center particle, and are centered around each particle; the black lines are constructed 
by finding the intersection of the green circles and indicate the width of the solid angle between the center particle and each of 
its neighbours, respectively. 



Combining Eqn. [T] and Eqn. [5] leads to a condition for 
the determination of the neighbor shell radius, 

Em 

- 7T~ < n,m+l, (o) 

m — I 

where R\ refers to the shell radius containing m parti- 
cles. To solve this inequality, we start with the smallest 
number of neighbors capable of satisfying Eqn.f2J m = 3, 
and increase m iteratively. During each iteration, we 
evaluate Eqn. [3] and the smallest m that satisfies the 
equation yields the number of neighbors Nf,(i) with i?- m ^ 
the corresponding neighbor shell radius. It is straight- 
forward to show that the algorithm converges, because 
the neighbor distance increases monotonically due to the 
sorting, and the cutoff radius de- 

creases monotonically, i?^ m+1 ^ < R^ . 

To highlight the differences and similarities between 
the geometry of the SANN algorithm and that of the 
Voronoi construction, we show in Figure[2]a 2d-schematic 
of both. The left panel shows the Voronoi construction 
that identifies all particles A to F as neighbors of the cen- 
ter particle. The Wigner-Seitz cell is indicated with black 
lines. Neighbor A shares only a small face and is fragile to 
thermal fluctuations. The right panel shows the SANN 
algorithm. To compare the Voronoi construction to the 
SANN algorithm, we make use of the fact that the algo- 
rithm is scale free, i.e. fij/R^ = (O.brij) / (O.bR^) 



where O.brij is simply the distance from particle i to the 

midpoint between i and j, and 0.5R\ is half the shell 
radius. The green circles have a radius equal to the half 
the shell radius of the center particle, and are centered 
around each particle; the black lines are constructed by 
finding the intersection of the green circles and indicate 
the width of the solid angle between the center parti- 
cle and each of its neighbours, respectively. Hence, one 
way to picture this method is to picture slowly grow- 
ing spheres around each particle. When the interescting 
planes associated with particle i (the black lines in the 
plot) yield solid angles summing to 4tt, the shell radius 
of particle i has been found. This proceedure is then 
repeated for each particle. In the schematic (Figure 2), 
Particle A is not a neighbor since there is no overlap be- 
tween the green circle around particle A and green circle 
around the center particle, hence the fragility problem 
highlighted in the discussion of the Voronoi construction 
is not present here. The black lines indicate the width of 
the solid angles and can be compared to the faces of the 
Wigner-Seitz cell. However, the faces are not identical to 
the real Wigner-Seitz cell. In general the faces are either 
larger or smaller than the real Wigner-Seitz faces. Note 
that, by definition, SANN extends each face to the shell 
circle (see Fig. [TJ, hence, the black lines in the SANN al- 
gorithm overlap sometimes, i.e. between particle D and 
E as well as E and F. 
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B. Algorithm 

Following the procedure outlined in the previous sec- 
tion we propose this simple scheme to determine the near- 
est neighbors of particle i: 

1. Compute distances n t j to all potential neighbors 
{j} from i. 

2. Sort possible neighbors {j} by their distance r^j in 
increasing order. 

3. Start with m = 3 (i.e. the minimum number of 
neighbors). 

4. Compute i? 4 M = £™ 1 r itj /(m - 2). 

5. If (ijf > r,. m+ i), then increment m by 1 and go 
back to step 4. 

6. Otherwise, m is Nb(i), i.e. the number of neigh- 
bors for particle i, and the associated neigh- 
bor shell radius. 

A C/Fortran implementation of the scheme can be 
found in the Supplementary Materials fl2l ]. 

C. Algorithm properties and extensions 

Before comparing our algorithm to the results from 
a fixed-distance cutoff and a Voronoi construction, we 
first discuss several inherent properties of our SANN al- 
gorithm and possible extensions to the algorithm. 

Pair-wise symmetry: For both a fixed-distance cutoff 
and the Voronoi construction, the neighbors are 
symmetric in the sense that if particle i is a neigh- 
bor of j, then j is also a neighbor of particle i. In 
SANN, this symmetry is not ensured, because every 
particle has its own neighbor shell radius. Thus the 
distance between both particles can be smaller than 
the shell radius of particle i and larger than that of 
particle j at the same time; hence, asymmetries can 
occur. However, we have found that the fraction of 
asymmetric neighbors is quite small: below 5% for 
the systems we studied. For many applications this 
might not matter, but in case it does there exists an 
easy way to make the algorithm " symmetric" : after 
having applied the standard algorithm to neighbor- 
ing particles i and j and having obtained R^™ 1 ^ and 
Rj , we compute rjj (the distance between i and 

j) and if both n t j < r\ and nj < Rj™\ then i 
is neighbor of j and j is neighbor of i, otherwise we 
remove the asymmetric neighbor from the set. 

Local volume: It is possible to assign a local volume 
to each particle. The Voronoi algorithm has as 
an obvious choice for the local volume the Wigner- 
Seitz cell, and by construction the sum of all local 



volumes adds up to the total system volume. For 
SANN one can think of many different definitions 
of a local volume, e.g. related to the shell or cutoff 
radius, but there is no inherent definition. Conse- 
quently, the sum of such local volumes does not by 
definition equal the system volume. Again, if it is 
important to attribute a volume to each particle, 
we can simply (and somewhat arbitrarily) rescale 
all volumes, such that their sum equals the total 
volume. 

Independence of space dimension: Although de- 
signed for three-dimensional space, we point out 
that the algorithm is valid without modification 
for any space-dimension with d > 2 (and for d = 1, 
it is obviously not needed) . In particular in higher- 
dimensional space, e.g. when studying the packing 
of hyper-spherical particles, easy-to-implement 
algorithms that go beyond the fixed-distance 
cutoff are scarce and SANN might be an attractive 
procedure. 

Next-nearest neighbors: In principle, the algorithm 
can easily be extended to yield a set of next-nearest 
neighbors, e.g. neighbor particles with a distance 
corresponding approximately to the second peak of 
the pair correlation function g(r). For this task the 
algorithm is performed twice as follows: in the first 
run, the nearest neighbors are computed without 
any modifications. Then all these nearest neighbors 
are discarded from the list of possible neighbors, 
and the algorithm is run a second time. Because 
the algorithm is scale-free, no modification to the 
algorithm is required, and the next-nearest neigh- 
bor shell is obtained. Note that simply increasing 
the total solid angle to 87r in Eqn.[5]does not work, 
as the solid angle contribution of the nearest neigh- 
bors would dominate due to the large shell radius 
R[ m \ As we shall see later in the paper, in prac- 
tice this extension does not work particularly well 
for finding next-nearest neighbors. 

Extension to a tunable algorithm The advantage of 
the current algorithm is that it is parameter-free. 
It assumes that the sum of all solid angles formed 
by neighboring particles adds up to 4-7T. However, 
exceptionally, there may be situations where this is 
choice is counter-intuitive - e.g. if a particle resides 
at a surface. In that case the solid angle 47r may 
be scaled by a factor a (although we do not recom- 
mend this practice). In that case, the cutoff radius 
R^' would be defined through: 

m m 

47ra = £>[l-ow(0 iJ )] = £>(1 -r^" ). (4) 

3=1 3=1 

The condition to stop at the current iteration then 
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becomes 

R { r ] = < n, m+1 . (5) 



III. SIMULATION DETAILS 

Below, we briefly describe the systems and the simu- 
lation methods used to produce the test configurations 
studied in this paper. Moreover, we provides details 
about the library used for the Voronoi construction and 
the implementation of our SANN algorithm. At the end 
of this section we briefly review the bond-order correla- 
tors that we use later to perform a structural analysis on 
some of the systems. 

In what follows, we denote the temperature by T, the 
pressure by P and the (number) density by p. The pack- 
ing fraction is defined as <p = ^pd 3 , where d is the par- 
ticle's diameter. In what follows, a will be the unit of 
length for both the hard-sphere and the Lennard-Jones 
systems whereas for carbon we will express the length 
in A. All distances presented in the manuscript will be 
expressed in the appropriate length units. 



A. Sample preparation 

Monodisperse hard-sphere configurations were pre- 
pared using an event-driven molecular dynamics simu- 
lation in an NVT ensemble with N = 86400 particles in 
a cubic box with periodic boundary conditions and with 
temperature T = 1, mass m = 1 and diameter d = 1. The 
system was equilibrated at a packing fraction (0 = 0.54) 
within the solid-liquid coexistence region (<frf = 0.492 
and 4> s = 0.543 [lJl) and at a higher packing fraction 
((f) = 0.61) beyond the hard-sphere glass-transition pack- 
ing fraction (4> g = 0.58). In both configurations only 1% 
of the particles are labeled as solid-like with a qe bond 
order criterion; note that the qe criterial will be discussed 
later in the text. For more details on these simulations, 
we refer the reader to Refs. [l4[ and (l5| . 

The carbon phases were simulated using the LCBOPI + 
potential [l6[ at the same conditions as the study on 
diamond nucleation in Ref. [l7| . namely P = 30GPa 
and T = 3750A for the 3-fold coordinated liquid and 
graphite, and P = SbGPa and T = 5000A' for the 4-fold 
coordinated liquid and diamond phases. Both conditions 
correspond to 25% under-cooling with a nucleation free- 
energy barrier equal to or larger than AG = 25/cbT (with 
ks Boltzmann's constant) preventing spontaneous crys- 
tallization of the metastablc liquid phase. All systems 
contained N — 1000 particles, with the exception of the 
graphite crystal which had N — 960 particles. More de- 
tails on the simulation methods and the semi-empirical 
interaction potential are given in Ref. (l6j . 

In our discussions of all Lennard-Jones systems, we 
denote with T* and P* the temperature and pressure 



in reduced units (T* = k B T/e and P* = Pa 3 /e), 
with e the Lennard-Jones well-depth, and p* = pa 3 the 
density in reduced units. To construct configurations 
of the Lennard-Jones fee crystal and liquid phases we 
performed Monte Carlo simulations in the isothermal- 
isobaric ensemble for particles interacting via a truncated 
and shifted Lennard-Jones pair potential [H, [l!| with 
a cutoff distance of 2.5. For both phases, a system of 
N = 4000 particles was prepared at the reduced tempera- 
ture T* = 0.92 and pressure P* = 5.68. Under these con- 
ditions, which correspond to 20% under-cooling with re- 
spect to coexistence, the liquid phase is metastable with 
respect to the fee crystal phase. However, a nucleation 
free-energy barrier of AG « 20fcsT prevents spontaneous 
crystallization on simulation time scales [2~oT]. The two- 
phase liquid-crystal system was simulated using the same 
Lennard-Jones potential at the same conditions, but with 
N = 8000 particles. The equilibration was biased with a 
quadratic potential on the number of solid-like particles 
to prevent the further growth of the crystal phase. See 
Ref. [l8| for details on biased Monte Carlo simulations 
and Ref. [1] on how to identify solid-like particles. The 
two-phase liquid- vapor configurations were prepared us- 
ing the same Lennard-Jones potential and equilibrated 
using NVT Monte Carlo simulations at reduced temper- 
ature T* = 1.0, number density p* = 0.3, and a system 
size of N = 5000 particles. In order to study the liquid- 
vapor interface, the simulation box was elongated along 
the x-axis such that the box length in the x direction was 
2.5 times as long as the y and z directions. As a result 
of this simulation box geometry, the resulting liquid-gas 
interface was perpendicular to the x-axis. 



B. Voronoi and SANN implementation details 

To compute the Voronoi construction, we used the 
Open Source Computational Geometry Algorithms Li- 
brary (CGAL [2l|), version 3.7. However, we were only 
interested in the set of nearest neighbors and not in the 
additional Voronoi information such as the volume, faces, 
edges, etc. of the Wigner-Seitz cell. Therefore, it was suf- 
ficient and computationally cheaper to perform a Delau- 
nay triangulation, which is the dual of the Voronoi con- 
struction. Either construction can be transformed into 
the other and both yield identical nearest-neighbor sets. 
We used CGAL's "exact predicates inexact construction" 
kernel and included 8 periodic copies of each particle to 
emulate periodic boundary conditions. Although CGAL 
does support 3d Dclaunay triangulation with 3d period- 
icity, it turned out that the run-time was significantly 
worse. The particles were inserted sequentially to map 
CGAL vertex handles and our particle ids. 

To speed up the SANN algorithm we made use of 
a Verlet list [l8| with a long cutoff distance to deter- 
mine the set of possible neighbors for each particle. Al- 
though this method involves a cutoff parameter, we could 
have chosen a parameter-free algorithm like a binary 
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space partitioning tree or an octree [221 ] . In general, any 
domain-decomposition method suffices as long as it pro- 
vides enough particles for the algorithm to converge. 



finish by presenting run-times of each algorithm for sev- 
eral simulation samples. 



C. Bond-order correlator 

In order to identify solid-like particles in some of the 
systems, we used local bond-order parameters according 
to Ref. Q. The orig inal order parameter described by 
Steinhardt et al. [23[ is based on the idea of expanding 
the neighborhood of each particle in a system in terms 
of a specific set of spherical harmonics, e.g. expanding 
in terms of the spherical harmonics with I = 3, I = 4 or 
/ = 6, depending on the local symmetry. The algorithm 
was later refined by ten Woldc et al. [J] for the study of 
nucleation, and has proven to be a useful tool even in the 
case of higher-dimensional systems [24|, [H[ . 

To compute the bond-order parameter each particle i 
is assigned a (21 + l)-dimensional complex vector qi(i) 
whose m-th component is defined by, 



(6) 



where Nb(i) denotes the number of nearest neighbors, 
Yim{fij) is the set of spherical harmonics of order I with 
components — I < m < I, f,j is the unit vector pointing 
from the center of i to its neighbor j, and the sum runs 
over all neighbors {j} of particle i. From this we can 
construct a measure for the neighborhood similarity of 
two particles, 



di(i,j) 



(7) 



where the superscript star denotes the complex conju- 
gate. We call the di(i,j) the local bond-order correla- 
tor, which is one when both particles are in an identi- 
cally ordered environment. To distinguish reliably be- 
tween solid-like and liquid-like particles, particularly in 
an under-cooled liquid, additional steps are required to 
increase the contrast. However, since a change in the 
neighborhood algorithm already affects this stage of the 
analysis, we will not follow the procedure to the end, but 
instead compare the local bond-order correlators. 



IV. RESULTS 

In what follows we apply the proposed algorithm 
(SANN) to several simulation samples and compare the 
resulting set of nearest neighbors to the sets obtained 
from both the fixed-distance cutoff criterion and the 
Voronoi construction. Moreover, on some systems we 
perform a structural analysis using bond-order param- 
eters and discuss the impact different nearest-neighbor 
sets have on the bond-order correlator distributions. We 



Bulk phases 

To start, we compute the nearest-neighbor distribu- 
tion P(N n ) for the bulk phases described in Section llll Bl 
using three neighborhood algorithms. For the fixed- 
distance cutoff, we set the cutoff to the minimum of 
the pair correlation function, which yields r c = 1.5 and 
r c = 1.35 for the Lennard- Jones liquid and fee crystal 
phases, respectively, r c = 2.0 for all carbon phases, and 
r c = 1.35 and r c = 1.3 for the low- and high-density 
hard-sphere suspensions. 

Figures [3^ and [5J3 depict for the liquid and fee 
Lennard-Jones system the nearest-neighbor distribution 
P(N n ) computed using fixed-distance cutoff distance 
(C), Voronoi construction (V) and our SANN algorithm 
for nearest neighbors (S) and for next-nearest neighbors 
(52). In both systems, the Voronoi construction identi- 
fies more nearest neighbors on average than the fixed- 
distance cutoff and SANN methods. Its peak in the 
nearest-neighbor distribution is around 14 neighbors for 
both the metastable-liquid and the fee phases. The fixed- 
distance cutoff (C) exhibits a nearest-neighbor distribu- 
tion which peaks around 13 neighbors in the liquid and 
at 12 in the fee crystal, whereas the distribution obtained 
using the SANN algorithm peaks around 11 — 12 neigh- 
bors for the liquid and sharply at 12 neighbors in the fee 
crystal. Note that 12 is also the number that one would 
expect from a close-packed arrangement of spheres. 

In order to get a better understanding regarding the 
particles which arc identified as nearest neighbors in the 
SANN and Voronoi algorithms, we also compute the pair 
correlation function using only the nearest neighbors ob- 
tained with each method (g n n(f)) and compare them to 
the g(r) of all particles. In the following discussions, we 
picture the environment around each particle to consist 
of several shells, each shell associated with a peak in the 
g(r). Hence, everything up to the first minimum of the 
g(r) corresponds to the first shell, everything between the 
first and the second minimum to the second shell and so 
forth. The fixed-distance cutoff method is not applied 
here as, by definition, its g n n(r) yields g(r) exactly up 
to the cutoff radius, after which it is zero. The upper 
graphs show the pair correlation functions g n n(r) and 
g(r) for reference, and the lower graphs show the ratio 
9nn{r) / g(r) . At a given distance r, the latter ratio gives 
1 if all particles at this distance arc identified as nearest 
neighbors, and reduces to zero if none of these particles 
are considered neighbors. Hence, a steep decrease in the 
ratio gnn{r) / g(r) indicates few fluctuations in the selec- 
tion of the neighbors. In addition to the nearest neigh- 
bors, the graphs also show results for the next-nearest 
neighbors obtained from the SANN method ($2)- 

Figure[3j: and[3ji plots these functions for the Lennard- 
Jones phases. They show that g n n(f) for the Voronoi 



7 





b) 






4.0 




3.0 


a 


2.0 


to 






1.0 




0.0 




1.0 


■£? 






0.5 


s 




to 


0.0 






v ■ 


.1 \ 


! s 

! s 2 ■ 




\ 

i i 


1 1.5 


2 2.5 



c) 



d) 



Figure 3: Nearest neighbors distribution P(N n ) for a Lennard- Jones liquid (panel a) and fee crystal (panel b) obtained by 
fixed- distance cutoff (C), Voronoi construction (V) and SANN considering neighbors belonging to the first coordination shell 
(S). Panels c) and d) plot the pair correlation functions g(r), considering all particles, as a reference (thin grey dotted line), 
and g n n(r), considering only nearest neighbors (Voronoi (V) and SANN (S)) and next-nearest neighbors (SANN (52)), for 
both the liquid (c) and the fee crystal (d). In addition, their fraction gnn(r)/g(r) is shown. 



construction (V) is identical to the reference g(r) up to 
the first minimum, and in the fee crystal even slightly 
beyond that. From the position of the decrease in the 
gnn(r)/g(r), i.e. slightly to the right of the first mini- 
mum in the g(r), we see that the Voronoi algorithm also 
includes some particles from the second neighbor shell 
(see upper panels of Figs. [3t an d EJi). This behavior 
originates from fluctuations which cause the Voronoi cell 
of next-nearest neighbors to occasionally share a small 
face @. There exist extensions to the Voronoi con- 
struction which attempt to increase the robustness of 
the algorithm to fluctuations @-[T3|. However, many of 
them introduce non-inherent parameters and, as such, 
are not parameter-free. Therefore, we will not consider 
them here. In contrast to the Voronoi algorithm, the 
gnnir) / g(r) associated with the SANN algorithm (S) 
drops to zero at the first minimum for both the liquid and 
crystal phases and therefore hardly includes any next- 
nearest neighbors. The SANN algorithm to determine 
next-nearest neighbors, denoted 52, does not yield very 
precise results. In particular, it identifies a fraction of 
particles in both the second and third shell of the liquid, 



and in the case of the solid even goes beyond the third 
shell. In both cases this can be attributed to the form 
of the g(r), i.e. the broadness of the second peak in the 
liquid, and the closeness of the second and third peaks in 
the solid. Unfortunately, this is a recurrent problem with 
trying to use SANN to determine next nearest neighbors, 
and for this reason 5*2 will not be discussed further in this 
paper. 

For both hard-sphere systems studied, i.e. (f> = 0.54 
(fluid) and <f> = 0.61 (glass), the nearest-neighbor distri- 
butions of Figures |3Ji andSJa computed using the Voronoi 
construction (V) present a peak around 14 neighbors, 
given that some of the neighbors from the second shell 
are included (as shown in Figs. 0fc and|3Ji). In contrast, 
the distribution obtained using the fixed-distance cutoff 
(C) algorithm and the one obtained with SANN (S) are 
fairly similar (Figs. H^, and[4}D) and both peaked around 
12 neighbors for both packing fractions. Again, this is the 
number one would expect in a close-packed arrangement 
of spherical particles. From FiguresSJ; andHJl we see that 
the pair correlation function for the Voronoi construction 
is identical to the reference g(r) up to the first minimum. 
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Figure 4: Nearest-neighbor distribution and g(r) as in Figure[3] but for a monodisperse hard spheres system at i 
a and c) and <j> = 0.61 (panels b and d). 
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But, as in the Lennard Jones system, it also seems to 
partially include particles from the second neighbor shell 
(see upper panels of Figs. 0t and d). In contrast, the 
gnn(f) computed using SANN drops to zero at the first 
minimum at both </>'s and does not include next-nearest 
neighbors. 

Figures [5] and [5] depict results for systems consisting of 
3-fold coordinated liquid carbon and graphite and 4-fold 
coordinated liquid carbon and diamond, respectively. In 
contrast to the Lennard- Jones and hard-sphere systems, 
carbon is a highly structured network-forming system, 
even in the liquid phase. It features open structures with 
few (3 or 4) close-by ordered neighbors: the 3-fold co- 
ordinated liquid carbon has a graphite-like structure in 
the first coordination shell, whereas the 4-fold coordi- 
nated liquid carbon has a rather pronounced diamond- 
like structure in the first coordination shell, shown in 
the strongly anisotropic angular distribution of the first 
neighbors. This is reflected in the pair correlation func- 
tions g(r) of Figures [SJ; and [Bfc, that show a sharp first 
peak followed by a broad deep minimum: a sign that up 
to the second neighbor shell the liquid has a structure 
almost as pronounced as the one of the corresponding 
solid. 

The upper panels of Figures [5] and [5] represent the 



P(N n ) computed for the 3-fold coordinated and 4-fold 
coordinated carbon phases, respectively. The nearest- 
neighbor distributions from the three methods differ sig- 
nificantly: with the cutoff distance set to the first min- 
imum of the <?(r), the fixed-distance cutoff yields a dis- 
tribution peaked sharply around 3 (Fig. [5j 3-fold coor- 
dinated system) and 4 (Fig. [51 4-fold coordinated sys- 
tem) particles both for the liquid and the crystal phases. 
Counting particles that form chemical bonds, the 3-fold 
coordinated carbon should have 9 neighbors within the 
first two shells (3 in the first shell, separated by a dis- 
tance of about 1.4 A, and 6 in the second shell belonging 
to the same graphite layer), whereas the 4- fold coordi- 
nated carbon should have 16 neighbors within the first 
two shells (4 in the first shell, separated by a distance of 
about 1.54 A, and 12 in the second shell). The SANN 
algorithm peaks around 10 (3-fold coordinated systems, 
Fig. [SJ) and 12 (4-fold coordinated systems, Fig. [5J par- 
ticles for both liquid and crystal phases. Those numbers 
indicate that SANN includes in each particle's neighbor 
list most of the particles belonging to the second coordi- 
nation shell, as confirmed by the g(r) plots in the lower 
panels of Figures [5] and [51 The Voronoi construction 
peaks around 17 (3-fold coordinated systems, Fig.[5j and 
20 (4-fold coordinated systems, Fig. [6j since it includes 
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Figure 5: Nearest-neighbor distribution and g(r) as in Figure [3l but for both a 3-fold coordinated carbon liquid (panels a and 
c) and graphite crystal (panels b and d). 



particles from the second and third shells. This happens 
in the liquids, graphite and diamond, as shown by the 
non-monotonic decay of the Voronoi's g n n(r)/g(r) that 
extends well beyond the second minimum in these cases. 

To explain the behavior for both the Voronoi construc- 
tion and the SANN method, we recall that both 3-fold 
and the 4-fold coordinated carbon phases are network- 
forming open structures. Moreover, graphite forms layers 
that are several particle diameters apart (with 3.4 A the 
distance between two layers). This structure affects both 
algorithms differently; by definition the Voronoi con- 
struction searches for neighbors that surround the cen- 
ter particle in all space dimensions, attempting to con- 
struct a complete (3d) Wigner-Setz cell. For the 4-fold 
coordinated carbon phases the Wigner-Seitz cell of the 
neighbors from the first coordination shell is a fragile 
tetrahedron, meaning that it is very likely that much 
further apart particles share a small face. Hence it con- 
tains particles from the second and even higher coordi- 
nation shells. The planar arrangement of neighbors in 
the 3-fold coordinated carbon phases forces the Voronoi 
construction to consider particles from neighboring lay- 
ers to complete a 3d Wigner-Seitz cell; particles that, as 
one might argue, belong to an entirely different neighbor- 
hood. In contrast, SANN does not attempt to complete 



a 3d environment. However, in order to complete its 
neighborhood with only neighbors from the first coordi- 
nation shell they must have almost identical distances. 
This is rarely the case in physical systems and therefore 
it includes more distant particles as well. Then, however, 
its neighborhood is dominated by the particles from the 
first coordination shell, since they are much closer and 
consequently contribute much larger solid angles in com- 
parison to more distant particles from the second shell. 
This effect seems to be more pronounced in the diamond 
crystal phase than in the other carbon phases, and it is 
the reason why not all particles from the second shell are 
included. We note that this effect can be mitigated by 
using the extension to a tunable algorithm as detailed in 
Section [II C[ or by giving neighbours a "weight" propor- 
tional to their solid angle. 

As an example for this behavior, we depict a few 
graphite layers in Figure [7J where we color particles iden- 
tified as neighbors with the Voronoi (panel a) and SANN 
(panel b) algorithms. Although geometrically correct, 
the additional neighbors from arguably different neigh- 
borhoods identified by the Voronoi construction may dis- 
tort results for local quantities. Figure [Ji; presents a top- 
view of the center layer of panel (b) and shows that the 
neighborhood identified by SANN includes the complete 
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first and almost complete second neighbor shell. 

Interfaces 

We now apply the three algorithms to two-phase sys- 
tems with planar interfaces, namely a liquid-crystal and 
a liquid-vapor Lennard- Jones systems as described in 
Section IIIIBI The two phases are arranged in a slab- 
geometry such that the interfaces are normal to the in- 
direction. For the fixed-distance cutoff we use r c = 1.5, 
which corresponds to the minimum of the g(r) for the 
liquid in the liquid-crystal system; note that the density 
of this liquid is not the same as the density of the liquid 
in the liquid-gas system. This choice of r c is arbitrary 
since there is no way to choose a cutoff which satisfies all 
four phases simultaneously. 

In Figure [5] we plot the average number of near- 
est neighbors (NNb(x)\ an d the corresponding variance 

Var( X ) = (N* m ( X )) - (N Nb (*)) 2 as a function of x for 
both systems. As expected, the liquid-crystal interface 
(Fig. [5] panel a) is barely visible from all three methods 
and the results are quite similar. However, while the 
fixed-distance cutoff and SANN methods seem to show 



a slight increase in {^N;sib{x)j in the crystal phase, this 
does not appear evident in the Voronoi algorithm. In all 
cases, the crystal seems to have a slightly lower variance 
than the fluid. In contrast, the liquid-vapor interface 
(Fig.[5]panel b) is very well captured by the average num- 
ber of nearest neighbors (upper panels) computed by the 
fixed-distance cutoff, whereas the SANN algorithm shows 
only a slight decrease of the number of nearest neighbors 
in the vapor phase. The Voronoi algorithm finds even 
more neighbors in the low-density vapor phase and its 
standard deviation (lower panel) increases strongly in the 
vapor phase. This behavior reflects the strong sensitiv- 
ity of the Voronoi construction to thermal fluctuations. 
Although Var(x) also fluctuates in the other two meth- 
ods, the changes are much less pronounced than in the 
Voronoi case. 

To get a better understanding for the nearest neigh- 
bors found at the interface and in the vapor phase, in 
Figure we show a snapshot from the two-phase liquid- 
vapor system where the vapor phase has been shifted to 
lie in the center of the box. We have selected three par- 
ticles (two at the liquid-vapor interface and one in the 
vapor phase) and have calculated their nearest neighbors 
using all three algorithms. As expected, in or at the va- 
por phase the fixed-distance cutoff (upper panel) finds 
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Figure 7: Simulation snapshot of 3-fold coordinated carbon graphite showing the first neighbors (yellow) of a center particle 
(gray). Surrounding particles that are not part of the neighborhood are shown in blue, a) Voronoi construction; b) SANN 
algorithm; c) Top-view on the center layer of panel b. 
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Figure 8: Results for two-phase samples in a slab-geometry, with the two interfaces oriented normal to the x-direction. The 
two phases are a) liquid-crystal and b) liquid-vapor. Note that the densities of the liquid in both cases are not the same. For 

both samples the upper panel shows, as function of x-position, the average number of nearest neighbors, ^iVjvf,(x)^ , and the 

lower panel the corresponding variance, Var(x) = (^Nf /b (x)^ — (^NMb(x)^ , for each of the algorithms: fixed-distance cutoff 

(C) with r c = 1.5, Voronoi construction (V) and SANN considering neighbors belonging to the first coordination shell (5). 
Note the different scales on the y-axis. 



few neighbors when the cutoff is not tuned to the va- 
por phase. Note that tuning the fixed-distance cutoff for 
the liquid and vapor phases simultaneously is not possi- 
ble. The Voronoi algorithm (center panel) detects many 
neighbors, both for particles at the interface and partic- 
ularly for particles in the vapor phase. At the interface 
the SANN algorithm (lower panel) finds neighbors mostly 
from the interface and liquid and, unlike the Voronoi con- 
struction, not far-off in the vapor. 



Application to bond-order parameters 

Finally, we apply the neighborhood algorithms to in- 
vestigate their effect on the local bond-order parameters 
used when studying crystal nucleation. To choose the 
order of the spherical harmonics in Eqn. [6j we match the 
symmetry of the spherical harmonics, i.e. I, to the sym- 



metry of the crystal under study. For the Lennard- Jones 
system we use I = 6 due to the close-packed crystalline 
structure of the fee (this is also what is typically used 
to study hard-sphere systems). In the original article on 
diamond nucleation [17| . I = 3 was applied to grow both 
carbon crystal phases since this order parameter is not 
able to distinguish between graphite and diamond struc- 
tures. The symmetry I = 3 was required since only the 
first neighbor shell was taken into account. However, in 
the present analysis, both Voronoi and SANN algorithms 
resulted in neighbor lists which included more than the 
first neighbor shell. Consequently, the symmetry of the 
neighborhood changes, and I = 6 becomes perfectly com- 
mensurate with the symmetry of this extended environ- 
ment. Therefore, we settled for I = 6 for all systems and 
set the fixed-distance cutoff in the carbon case to 2.7 (the 
minimum after the second peak of the g(r)) to include 
next-nearest neighbors. The results for the local bond- 
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Figure 9: Visual representation of the three algorithms in a 
two-phase liquid-vapor Lennard- Jones system: fixed-distance 
cutoff (top), Voronoi construction (center) and SANN (bot- 
tom). In each case, we select the same particles and check 
which neighbors are detected using each algorithm. 



order correlators distribution P[de(i,j)] are presented in 
Figure 021 

The main criterion for a good order parameter is to 
have as little overlap as possible between the crystal 
phase (upper panels) and the (meta-stable) liquid phase 
(lower panels). As shown in Figure IT0"k all neighborhood 
algorithms perform reasonably well for the Lcnnard- 
Jones system. However, for the graphite crystal (Fig.llOb 
upper panel) both SANN and the fixed-cutoff algorithms 
allow one to distinguish between a liquid and a solid 
environment, whereas the Voronoi construction shows 
a severe broadening of the distribution of the graphite, 
causing a substantial overlap with the liquid's distribu- 
tion (lower panel) . With such an overlap the bond-order 
parameters would fail to distinguish between liquid-like 
and solid-like particles. In the diamond case (Fig. [TOb ) 
the SANN algorithm performs worse than the others 
and again the fixed-distance cutoff works best. In the 
graphite case the failure of the Voronoi construction can 
be attributed to the (arguably) spurious neighboring par- 
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Figure 10: Distribution of the local bond-order correlator 
d§(i,j) using different neighbor criteria. Panel a) shows re- 
sults for Lennard- Jones fee crystal (upper panel) and liquid 
phases (lower panel), panel b) for the 3-fold coordinated car- 
bon graphite (upper panel) and liquid (lower panel) phases, 
and panel c) for the 4-fold coordinated carbon diamond (up- 
per panel) and liquid (lower panel) phases. For the carbon 
phases the fixed-distance cutoff distance was set to 2.7 to in- 
clude the next-nearest neighbors, as do inherently both the 
Voronoi and SANN algorithms. 
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System 


C 


V 


S 


v/c 


s/c 


Lennard- Jones liquid 


25 


610 


45 


24.4 


1.8 


Lennard- Jones fee 


20 


753 


48 


37.7 


2.4 


Hard-Spheres 4> = 0.54 


507 


14390 


1022 


28.4 


2.0 


Hard-Spheres <p = 0.61 


528 


15050 


1091 


28.5 


2.1 


Carbon 3-fold liquid 


5 


160 


8 


32.0 


1.6 


Carbon 3-fold graphite 


5 


190 


7 


38.0 


1.4 


Carbon 4-fold liquid 


6 


153 


10 


25.5 


1.7 


Carbon 4-fold diamond 


6 


180 


9 


30.0 


1.5 



Table I: Run-times in milli-seconds of the fixed-distance cutoff 
(C), the Voronoi construction (V) and the SANN algorithm 
(S), and their ratios V/C and S/C. For implementation details 
we refer to Section IIII Bl for details on the system samples 
to Section fill Al and for the benchmarking procedure to the 
main text. 

tides located in different graphite layers, as previously 
discussed (Fig. [7]) . The behavior of the SANN algorithm 
in the diamond case finds its origin in that not enough 
neighbors from the second neighbor shell are included 
(on average a total of 12 instead of 16), and bond-order 
parameters are particularly sensitive to missing or ad- 
ditional neighboring particles. As already pointed out 
before, this effect can be mitigated by using the exten- 
sion to a tunable algorithm as detailed in Section III CI 
For example, using an a = 1.5 in the diamond case re- 
sults in the inclusion of all neighbor up to the complete 
second coordination shell (i.e. 16) and a P[de(i,j)} dis- 
tribution almost identical to that of the fixed-distance 
cutoff (data not shown). Alternatively, as mentioned be- 
fore, one could assign neighbours a "weight" proportional 
to their solid angle. 

Benchmarking 

Finally, we measure the run-time of each algorithm 
for all bulk phases discussed previously The benchmark 
was performed on a computer equipped with an Intel 
Core2 Quad Q9550 processor running at 2.83 GHz and 
4 GB of DDR2 RAM running at 1066 MHz. All source 
code was compiled using the GNU gec compiler version 
4.5.1. The operating system was a 64-bit OpcnSuSE 
Linux with Kernel 2.6.37. Each algorithm was running 
single-threaded and computed the neighbor sets for all 
particles in the system. To improve accuracy whenever 
the run-time was near our time resolution we measured 
the total time of 10 sequential repetitions. All data pre- 
sented are averages over at least 3 separate program runs. 

The timings are presented in Table HI Compared to the 
fixed-distance cutoff the Voronoi construction takes 24.4 
to 38.0 times longer to compute. In contrast, the com- 
putational cost of SANN is only 1.4 to 2.4 times that 
of the fixed-distance cutoff and thereby outperforms the 
Voronoi construction by an order of magnitude. There- 
fore, we consider SANN well-suited for application on- 
the-fly in simulations. 

As a final remark we like to point out that timing re- 



sults are highly implementation dependent and as such 
should be considered as indications only. On the one 
hand, the CGAL library used for the Voronoi construc- 
tion is reasonably fast, but faster implementations may 
be available. On the other hand, our implementations of 
the fixed-distance cutoff and the SANN algorithm may 
have room for optimization, too. 



V. CONCLUSION 

In this paper we have described an algorithm to com- 
pute a particle's nearest neighbors in an arbitrary many- 
particle system. The algorithm is similar to a fixed- 
distance cutoff in that all particles within a cutoff dis- 
tance are considered nearest neighbors. But rather than 
using one cutoff for all particles, this algorithm assigns to 
each particle an individual cutoff distance, thereby mak- 
ing it suitable for systems with inhomogencous densities, 
such as gravitational or multi-phase systems. The cutoff 
distance follows from a geometric requirement, namely 
that the sum of all solid angles associated with neighbor- 
ing particles adds up to Att. Thus, the algorithm becomes 
parameter-free and scale-free. Though the approach was 
inspired by the Voronoi construction, it has several ad- 
vantages over it: the presented algorithm is significantly 
easier to implement, computationally less expensive and 
more robust against thermal fluctuations. 

We tested the algorithm on a number of bulk phases in- 
cluding supercooled liquid and crystal phases of Lennard- 
Jones particles, hard spheres and 3-fold and 4-fold coordi- 
nated carbon. We compared the nearest-neighbor distri- 
butions obtained from SANN to both the fixed-distance 
cutoff criterion and the Voronoi construction. In the case 
of the Lennard- Jones and hard-sphere phases, our algo- 
rithm reproduces very well the nearest-neighbor distri- 
bution of a well-tuned fixed-distance cutoff. This is in 
contrast to a Voronoi construction, which has large fluc- 
tuations and as such does not perform as well. For the 
carbon phases, our algorithm includes the second neigh- 
bor shell, like the Voronoi construction, but avoids neigh- 
bors in the neighboring graphite layers. 

We also examined particles at the interface of two- 
phase systems, such as Lennard-Jones liquid-vapor and 
liquid-crystal systems. We find that when two high- 
density phases coexist, all algorithms give reliable re- 
sults. However, at the interface between a fluid and a 
low-density vapor, our algorithm is more robust to ther- 
mal fluctuations than the Voronoi construction. 

We then employed the neighbor information of all al- 
gorithms as input for a bond-order analysis, which is typ- 
ically used in crystal nucleation studies for the identifi- 
cation of solid-like particles in a supercooled metastablc 
liquid. Comparing the bond-order correlator distribu- 
tions, we found little difference between the algorithms 
for the Lennard-Jones system, indicating that all algo- 
rithms are suitable for structure analysis of close-packed 
systems. However, the Voronoi construction failed for 
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the graphite phase due to the identification of spurious 
neighbors located in different graphite layers, and SANN 
performed poorly in the diamond phase, due to the fact 
that not enough neighbors from the second neighbor shell 
were included. Hence, care has to be taken when apply- 
ing either one SANN or the Voronoi construction to open 
structures and network formers. But where the Voronoi 
construction fails SANN might succeed, and vice-versa. 

Finally, we performed benchmarks on the run-time for 
all algorithms. On all systems tested we found the com- 
putational cost of SANN to be at most 2.4 times that of 
the fixed-distance cutoff and in all cases it outperformed 
the Voronoi construction by at least an order of magni- 
tude. 

To conclude, when studying a system at several con- 
centrations or a heterogeneous system, the proposed al- 
gorithm has the advantage that is does not require tuning 
a parameter for every concentration/environment. Given 
the robustness and low computational cost of our algo- 
rithm, we argue that SANN is well suited not only for 
post-analysis, but also on-the-ffy in simulations. It reli- 
ably identifies the nearest neighbors, and its behavior for 
graphite and at a two-phase interface suggests its appli- 
cation to situations where the Voronoi construction suf- 
fers from distorted polyhedra, like in structural analysis 



of protein folding trajectori es [26l |27| , in DNA- mediated 
colloidal crystallization [2^, l29l|. in suspensions of patchy 
colloids with tetrahedral or octahedral symmetry [30l - l33j 
and in water [34| . Finally, the SANN algorithm is not 
only useful for simulation data, but should also be use- 
ful in analysing experimental 3D images, as obtained, for 
instance, by confocal microscopy or by tomography. 
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